module setlammod
	use maxRPmod
	use compute
	
	implicit none

	real(4), allocatable		:: LR(:,:,:)
	real, allocatable			:: xlist0(:,:)
	
	contains

	subroutine setlam0
		integer, parameter	:: nbase=10
		integer	:: j
		real	:: x(nx0)
		
		n0active=0
		if(allocated(lam)) deallocate(lam)
		if(allocated(LR)) deallocate(LR,xlist0)
		allocate(lam(0))
		allocate(LR(0,0:1,nsimeff),xlist0(nx0,n0max))
		call savelam
	end subroutine 
	
	function thdist(th1,th2) result(val)
		real	:: th1(nth),th2(nth),val
		val=min(norm2(th1-th2),norm2([th1(4:6),th1(1:3)]-th2))
	end function
	
	subroutine setlam(idirty)
		integer				:: idirty
		integer, parameter	:: nbase=1, ntries=50,ncollect=15, nfinal=200
		integer, parameter	:: niter0=50
		integer				:: n0new, nprevious
		integer				:: j,l
		real				:: RP, th(nth), x(nx0)
		integer				:: ntry, niter
		real,allocatable	:: cRPs(:),cxs(:,:),cths(:,:)

		print *,"XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
		print *,"XXXXXXXXX  entering setlam   XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
		print *,"XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
		niter=niter0
fpoint:	do 		
			ntry=0; n0new=0; nprevious=0
			cRPs=[real::]
			if(allocated(cths)) deallocate(cths)
			allocate(cths(nth,0))
			do
				ntry=ntry+1
				call rnun(x)
				call settailflag

				call maxfromx(niter,x,RP)
				call setth0(x,th)
print *,tailflag,RP,getRPstddevclust(th)				
				if(RP>level) then
					j=0
					do 
						j=j+1
						if(j>n0new) exit
						if(thdist(cths(:,j),th)<0.1) exit
					enddo
					if(j>n0new) then	
						n0new=n0new+1
                        thlist0(:,n0active+n0new)=th
						xlist0(:,n0active+n0new)=x
						cRPs=[cRPs,RP]
						cths=reshape([cths,th],[nth,n0new])
						print *,"new point for RP=",RP
						call mdisp([x,th])
					else
						print *,"point previously identified"
                        nprevious=nprevious+1
					endif
				endif
				if(n0new>=ncollect .or. (n0new>0 .and. ntry>ntries) .or. (nprevious>3) ) exit
				if(n0new==0.and.ntry>nfinal ) exit fpoint
				if(idirty>0 .and. ntry>20 .and. maxval(cRPs)<level) exit fpoint
				if(ntry>ntries .and. niter==niter0) then
					niter=100
					print *,"switching to 100 iterations searches"
					cycle fpoint
				endif
			enddo
			print *," after ",ntry," searches, adding ",n0new, " points for a total of", n0active+n0new
			call mdisp(cRPS)
			call addpoints(n0new)
		enddo fpoint
	end subroutine


	subroutine addpoints(n0newo)
		integer				:: n0newo, n0new
		integer				:: l,j,lb
		integer				:: l1,lb1,l2,lb2
		real(4), allocatable	:: LRnew(:,:,:)
		real, allocatable	:: f0s(:,:,:,:)
		real				:: th(6)

		if(lowramflag==.false.) then
			n0new=n0newo
			allocate(LRnew(n0active+n0new,0:1,nsimeff))
!$omp parallel do num_threads(10)
			do l=1,nsimeff
				LRnew(1:n0active,:,l)=LR(:,:,l)
			enddo
			call move_alloc(LRnew,LR)
		else
			n0new=n0active+n0newo
			n0active=0	
			deallocate(LR)
			allocate(LR(n0new,0:1,nsimeff))
		endif
		
		allocate(f0s(2,n0new,npb,nblocks))
!$omp parallel do private(l,j,th) collapse(2)		
		do lb=1,nblocks
			do l=1,npb
				do j=1,n0new
					th=thlist0(:,n0active+j)
					f0s(1,j,l,lb)=getdens0(th(1:3),Y0s(:,l,lb))/densimp(l,lb)
					f0s(2,j,l,lb)=getdens0(th(4:6),Y0s(:,l,lb))/densimp(l,lb)
				enddo
			enddo
		enddo
		if(stage==2) then
!$omp parallel do private(l1,lb1,l2,lb2) num_threads(10)
			do l=1,nsimeff
				l1=pairsind(1,l)
				lb1=pairsind(2,l)
				l2=pairsind(3,l)
				lb2=pairsind(4,l)
				do j=1,n0new
					LR(n0active+j,0,l)=0.5*(f0s(1,j,l1,lb1)*f0s(2,j,l2,lb2)+f0s(2,j,l1,lb1)*f0s(1,j,l2,lb2))
				enddo
			enddo
		else
!$omp parallel do private(l1,lb1,l2,lb2) num_threads(10)
			do l=1,nsimeff
				l1=pairsind(1,l)
				lb1=pairsind(2,l)
				l2=pairsind(3,l)
				lb2=pairsind(4,l)
				do j=1,n0new
					LR(n0active+j,0,l)=f0s(1,j,l1,lb1)*f0s(2,j,l2,lb2)
				enddo
			enddo
		endif
			
!$omp parallel do private(j)
		do l=1,nsimeff
			do j=n0active+1,n0active+n0new
				LR(j,1,l)=getdens0p(thlist0(:,j),pairsdat(l))
			enddo
		enddo
		n0active=n0active+n0new
		lam=[lam,[(0.05,j=1,n0newo)]]
		call adjustlam
	end subroutine
	
	subroutine adjustlam
		real, dimension(n0active)		:: lamfac, RPold,RP
		real(4), dimension(n0active)	:: mlam,RPc
		integer		:: lc,i,l,j
		integer, allocatable		:: isel(:)
		real, allocatable			:: lamdum(:)
		real						:: leveltarget
		real(4), allocatable		:: LRnew(:,:,:)
		real, allocatable			:: infostats(:,:)
		
		leveltarget=level-levelfudge
		lamfac=4
		RPold=leveltarget
		RP=leveltarget
		do lc=1,5000
			if(lc==300) print *,"reached 300 iterations"
			if(mod(lc-2,50)==0) then
				call mdisp([real(lc),maxval(RP)])
				print *,"RP and lam"
				call mdisp(RP.cvr.lam)
			endif

			if(((lc>200).and.maxval(RP)<leveltarget+0.4*levelfudge) .or. (lc>500)) then
				exit
			endif
			if(lc<400) then
				lam=lam*exp(lamfac*(RP-leveltarget))
			else
				where(RP>leveltarget+0.1*levelfudge)	 lam=lam*exp(lamfac*(RP-leveltarget))
			endif
			where((RP-leveltarget)*(RPold-leveltarget) >= 0.0)
				lamfac=lamfac*1.1
			elsewhere
				lamfac=lamfac*0.5
			endwhere
			
			call setbounds(lamfac,0.001,20.0)
			call setbounds(lam,1E-10,10000000.0)
			RPold=RP
			mlam=lam
			RPc=0
!$omp parallel do reduction(+:RPc)
			do l=1,nsimeff
				if(sum(LR(:,1,l)*mlam)<cvglobal) RPc=RPc+LR(:,0,l)
			enddo
			RP=RPc/(nsim*npb)
		enddo
		call mdisp(RP.cvr.real(lam))
		print *,minval(RP),maxval(RP)
		lam=merge(lam,0.0,.not. ( (lam<1E-4 .and. RP<leveltarget-0.4*levelfudge) .or. lam<1E-8))
		call setctest
		
		isel=pack([(i,i=1,n0active)],.not. ( (lam<1E-4 .and. RP<leveltarget-0.4*levelfudge) .or. lam<1E-8))
		
		if(size(isel)<n0active) then
			print *,"keeping ",size(isel), " points, deleting ",n0active-size(isel)
			lamdum=-lam(isel)
			call svrgp(lamdum,lamdum,isel)
			lam=lam(isel)
			n0active=size(isel)
			thlist0(:,1:n0active)=thlist0(:,isel)
			xlist0(:,1:n0active)=xlist0(:,isel)
			if(lowramflag==.false.) then
				allocate(LRnew(n0active,0:1,nsimeff))
!$omp parallel do num_threads(10)
				do l=1,nsimeff
					LRnew(:,:,l)=LR(isel,:,l)
				enddo
				call move_alloc(LRnew,LR)
			endif
			print *,"done thinning"
		endif
		call mdisp(xlist0(:,1:n0active))
		call mdisp(thlist0(:,1:n0active))

		call mdisp(RP(isel).cvr.real(lam))
		print *,minval(lam),maxval(lam)
		call savelam		
		call printtime
	end subroutine

	subroutine setctest
		integer	:: lb,l1,i2,ix,l
		integer(8)	:: ccode,tcode
		real(4)		:: mlam(n0active)

		mlam=lam
!$omp parallel do private(l,l1,i2,ix,tcode,ccode)
		do lb=1,nblocks
			l=pairbind(lb)
			do l1=1,npb
				do i2=1,nbi
					tcode=0
					ccode=credflags(i2,l1,lb)
					do ix=1,nfpc
						tcode=2*tcode
						if(ccode<0) then
							l=l+1
							if(sum(LR(:,1,l)*mlam)<cvglobal) tcode=tcode+1
						endif
						ccode=2*ccode
					enddo
					ctest(i2,l1,lb)=tcode*2**(64-nfpc)
				enddo
			enddo
		enddo
	end subroutine
end module